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to that presented by I. Tarnove (3), but with revisions in the root-finding and 
scaling procedures. The capability of calculating the associated eigenvectors 
has also been added to the previously-published algorithm. 

This method uses an iterative root-finding technique and will converge 
with almost any first estimate. To preserve accuracy, the original matrix is 
used throughout for the evaluation of the function and scaling is done by powers 
of two. 

A FORTRAN-63 program for the CDC-1604, using this method, has been written 
for matrices having polynomial elements. The code was evaluated using a number 
of trial problems, including pathological matrices, and some of these results 
have been incorporated into this report. 



Prepared by: R. D. Brunell 

Mathematician, Computer Facility 

Approved by: Released by: 

D. G. Williams C. E. Menneken 

Head, Computer Facility Dean of 

Research Administration 

U. S. Naval Postgraduate School Technical Report/Research Paper No. 69 

July 1966 



UNCLASSIFIED 



TABLE OF CONTENTS 



Title Page 1 

I. Statement of the Problem 2 

II. Method of Solution 3 

III. General Discussion 5 

IV. Numerical Examples 6 

V. Areas of Further Investigation 15 

References 16 

Appendix I: Flow Chart 17 

Appendix II: Subroutine Listings 20 

GENEIG 20 

CALNEWX 23 

CDTERM 25 

JORCOM 27 

CALFUNC 29 

CKEVEC 30 

FUNCEV 31 

ROMAT 32 

ACCEPT 33 

SCALE2 34 

SCALE 1 36 

EVER 37 

Appendix III: Subroutine Description 41 



F4-NPGS-GENEIG, A Method for the Solution of the 
"Generalized Eigenvalue Problem" with Polynomial 
Elements . June 1966. 



AN ITERATIVE SOLUTION TO THE GENERALIZED EIGENVALUE - EIGENVECTOR PROBLEM 



I. Statement of the Problem 

The "generalized eigenvalue problem" may be stated as follows: 

Find those values of the complex parameter, z, which will provide 
a non-trivial solution to the system of linear equations, 

A^^(z)X^ + + ... + A^^(z)X^ = 0 

A2^(z)X^ + A22(z)X2 + ... + ° 

• • • 

• • • 

A^,(z)X^ + A (z)X, + ... + A (z)X = 0. 
nl 1 nz z nn n 

This can be written more concisely as 
(2) A(z)X = 0 

where A is an n x n matrix, the elements of which are arbitrary complex 
functions of z, and X is a complex column vector of n elements. The set 
£z^:A(z^)X^ = o} is the set of eigenvalue^ and the X^ are called the 
associated eigenvectors. 

The special case where the A^^(z)’s are polynomials in z of degree 
2 or less is encountered most often in practical work. The rest of the 
discussion in this paper will be confined to problems of this type. However, 
it should be pointed out that the general method described here is not so 
restricted. 
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Under the above restriction we can rewrite A(z) as 



( 3 ) 



A(z) = A z^ + A-Z + A_ 
o i Z 



where the A^’s are constant complex matrices. We wilj. need this form 



in the later discussion on scaling. 

II . Method of Solution 

The heart of the method described here is the iterative root-finding 
technique. After a literature search, in which much background information 
was gained from Wilkinson (1), and some experimentation, it was decided that 
a method developed by J. F. Traub (2) would be used in this code. Muller’s 
method (7) has been mentioned by a number of authors in connection with 
iterative solutions in the complex plane; however, Traub ’s technique was 
chosen because of its simplicity and ease of computation. The methods are 
similar; both are of second order and have the same as 3 miptotic error constant. 
The main difference is that Traub uses the Newtonian form of the interpolating 
polynomial and Muller uses the Lagrangian form. 

In essence the technique is to perform a quadratic fit using the last 
three iterants and their functional values. This algorithm is defined as 
follows : 



(4) 




where 



(5) 




and 



( 6 ) 




The variable sign in equation (4) is chosen so that the denominator has 



the larger magnitude. 
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Convergence is determined by either of the two tests. 



( 7 ) 






or 



( 8 ) 




If either of these inequalities is met the process is considered to 
have converged to a root. 

In order to avoid convergence to a previously-calculated root, a 
method of suppressing zeros, first suggested by G, E. Forsythe (Reference 3, 
page 164) was incorporated into the code. The first root, is found 



One of the problems encountered in a program of this type is floating- 
point exponential overflow when working with matrices of large order. To 
alleviate this problem an automatic scaling feature was devised. This 
scaling is done in two parts, ^ priori and ^ posteriori . The ^ priori 
scaling is done at the beginning of the problem and is carried out in the 
following manner. The coefficients of each polynomial element in a row of 
the matrix A are summed and-a maximum is found. In a binary computer this 
maximum may be expressed exactly as bx2’^,. Then each coefficient in every 
polynomial in that row is divided by 2^ , This operation is carried out 
for each row in the matrix, 

A posteriori scaling is done only when an overflow occurs. If, at 

the time of the overflow, the estimate to the root being calculated is some 
k 

number, say cx2 , then A(z) (defined in equation (3) ) is replaced by 



as a zero of f(z) and the remaining roots, z^, are found as zeros of 



g^(z), where 



i 



(9) gj^(z) 5 f(z)/TT , i> 1 . 



k=l 
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( 10 ) 



A»(z») 



= A z 



o 




z + 



( 2 “)^ 



> 



and a root, z\ is then found for the new problem. The root of the 
original problem is then calculated as 
( 11 ) 2 . = 



The scaled problem is started by taking the last three iterants in the 

k 

original problem and dividing them by 2 . After a root is found the 
matrix is returned to its original form and the program proceeds to find 
the next root. 

In many cases it is desirable to obtain the eigenvectors associated with 
the eigenvalues. An option to perform these calculations is available 
in the program. Because there exists no procedure which will absolutely 
assure the finding of an eigenvector, the code will make 2n attempts, 
where n is the order of the matrix, at computing an eigenvector for each 
eigenvalue. This procedure is carried out as outlined below. 

An assumption is made that the elements of the vector may be written 
as multiples of a single element in the vector. One possible representation 
of this is 



( 12 ) 



T 

X. = (x b-x b .) 

1 li, 1 2i, 2 3i, n-1 ni 



Using this assumption the code makes its first attempt by setting 
and solving the system of equations (13): 



( 13 ) 



'^22 ^ 23 * • • ^ 2 n 



^32 ^33’ ' ' ^3n 



a - a . . .a 

n2 ni nn 



hi] 






hi 


= 


-^31 

• 


\hii 




• 

# 

\hil 
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where a., is the function A., evaluated at z.. 

ij 1 

The expression, A(z) X, is then calculated using this vector 
and the result compared to the zero vector. If the vector, X, is 
found to be unacceptable a second attempt is made by replacing a^j^ by a 
k=l,2,...,n, in (13). If this fails then is set equal to unity 
and the system (14) is solved, 



(14) 



/ 



ail ^13' ' • ^In 



\ 



^31 ^33 ’ • '^3n 



\ a - a ^ ... a / 
' nl n3 xml 

and checked as indicated above. 



1 X \ 




—a 


I 




12\ 


X 




— a_ 1 


3i 




32 1 


Vnil 




r^n2/ 



A fourth attempt may be made by replacing. 



in (14), a^j^ by for k=l,2,.. 
until an eigenvector is found. 



.,n. The process continues for all X.. or 

Ji 

If no eigenvector is found a statement to 



that effect is printed along with the results of the final attempt. 



Ill . General Discussion of Method 

A code, using the method described in this paper, has been written 
in FORTRAN-63 for use on a CDC 1604. This program has been checked out 
on a variety of problems, several of which are shown in the Section IV. 
The accuracy of the process has been shown to be quite good. This may 
be attributed to several factors. First, only the original input matrix 
is used for the evaluation of the function each time, and, second, all 
scaling is done by powers of two ( Reference (8) , page 98 ) . The method 
has been shown to be reasonably fast, although a good relative comparison 
is impossible because the number of iterations required is a function of 
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the estimate provided. Unfortunately this information generally has not 
been included in previous literature on this subject. Because the number 
of computations necessary to determine the next estimate to the root is 
so small compared to those required to evaluate the determinant, the time 
required for each iteration can be considered substantially equal to the 
time necessary for the latter. The speed of the program is further enhanced 
by the fact that through the use of good estimates only eigenvalues of 
interest need be computed, although convergence seems to result from almost 
arbitrary estimates. Also because of the feature which forces avoidance of 
previously found zeros, it is possible to restart the program and by-pass 
any such prior calculations. 

IV. Numerical Examples 

Most of the examples used for this section have been taken from papers 
on the standard eigenvalue problem. This was necessary because of the 
lack of comparative data on the generalized problem. However, these examples 
should be adequate to demonstrate the ability of this method to find path- 
ological roots and also to get some comparison as to its accuracy. All 
problems, using the method described in this report, were run on a CDC 1604 
which has a 48-bit word length. 
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Example 1 



(See Parlett (5) ). These results were obtained using Laguerre*s 



method on 


an IBM 7090 


which has 


a 36-bit 


word. A is 


defined as 






j 611-z 


196 


-192 


407 


-8 


-52 


-49 


29 \ 


196 


89 9- z 


113 


-192 


r-71 


-43 


-8 


-44 \ 


-192 


113 


899-2 


196 


" 61 


49 


8 


52 1 


407 


-192 


196 


611-2 


8 


44 


59 


-23 


-8 


-71 \ > 


61 


8 


411-2 


-599 


208 


208 


-52 


-43 


49 


44 


-599 


411-2 


208 


208 


-49 


-8 


8 


59 


208 


208 


99-2 


-911 


y 29 


-44 


52 


-23 


208 


208 


-911 


99-2j 






Table 


of Eigenvalues 
















. 









Actual (8 digits) 


Parlett (5) 


Present Method 


1020.0490 


1020.0500 


1020.0490 


1020.0000 


1019.9997 


1020.0000 


1019.9020 


1019.9019 


1019.9019 


1000.0000 


1000.0001 


1000.0000 


1000.0000 


999.99999 


1000.0000 


.098048640 


.098045509 


0.098048646 


0.0 


-0.0000094 


-0.000000003 


-1020.0490 


-1020.0490 


-1020.0490 



Total running time: 16 seconds 

Average number of iterations per eigenvalue: 10.5 
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Example 2 



(See Barlow and Jones (4) ) . Their work was done on an IBM 7040 

(36"bit word) using an extension of the classical secant method to the 
complex field. This paper (4) also gives results obtained by Eberlein (6) 
using a Jacobi-like method. 



A = 



11 


6 


-9 


-15 


3-z 


9 


—3 


-8 


6 


6-z 


-3 


-11 


7 


5 


-3-z 


-11 


12 


5 


-10 


-16- 



\ 



/ 



Table of Eigenvalues 



- Actual (6) 


Eberlein (6) 


Barlow & Jones (4) 


Present Method 


-Ko 

1.50016 + 3.570641 
1.50016 + 3.570641 
1.50016 - 3.570641 
1.50016 - 3.570641 


-.999 

1.505 + 3.571 
1.495 + 3.571 
1.505 - 3.571 
1.495 - 3.571 


-1.0000004 
1.49850 + 3.570271 
1.50150 + 3.571101 
1.49996 - 3.571271 
1.49864 - 3.570231 


-.999999999 
1.50001 + 3.570731 
1.50003 + 3.570691 
1.50001 - 3.570741 
1.49999 - 3.570701 



Total running time: 4.38 seconds 

Average number of iterations per eigenvalue: 6 
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Example 3 



(See Parlett (5) ). This example was used to demonstrate the 

ability of the code to extract only a few desired roots. 



-( 



c 

4C 



2C 

3C 



A = (B - zl) 







/-2 


2 


2 


2\ 


\ c -1^^ ~ 


d\ 




3 


2 


2 \ 


1 ^ ■ 1 5D 




0 


4 


2 


1 




\-l 


0 


0 


5/ 



The true eigenvalues of this system (Eberlein (6) ) are 60+20i, 

45+15i, 3p+10i, 15+5i, -12+4i, -9+3i, -6+2i, -3+i. 



Initial Guess 


Eigenvalue Found 


-2.5 - 0.625i 


-3.0000006 - 9.99999471 


-2.5 + 0.625i 


-3.0000005 + l.OOOOOOOi 


-5.0 - 1.25i 


-6.0000055 - 1.99999481 


-5.0 + 1.25i 


-5.9999975 + 2.00000031 



Total running time: 55.28 seconds 

Average number of iterations per eigenvalue; 17.5 
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Example 4 



In this example both the eigenvalues and eigenvectors were found. 



A = 




-2 



3 



-2 



\ 



1 



5-z -1 -1 




3 2-z -2 

-2 6 -3-z I 



Table of Eigenvalues and Eigenvectors 



True Values 


Computed Values 


Eigenvalues 


Eigenvectors 


Eigenvalues 


Eigenvectors 


-1 


1 


-1.00000000 


1.0 




0 




0.0 




0 




0.0 




1 




1.0 


2 


1 


1.999939 


1.0 




-6 




-5.9993 




-9 




-8.9991 




-8 




-7.9993 


2 


1 


2.00007 


1.0 




-6 




-6.0007 




-9 




-9.001 




-8 




-8.0008 


2 


1 


1.99998 


1.0 




-6 




-5.9998 




-9 




-8.9998 




-8 




-7.9998 



Total running time: 3.83 seconds 

Average number of iterations per eigenvalue: 8.25 
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Example 5 



This is an example of the generalized problem. 



2z +14Z-88 



-z -13Z-22 



z +2z-99 



A = 



-3z +21Z-36 



6z -6z-36 



4z -48z+108 



8z +8z-160 



z +7z+10 



-z +4z+45 



Table of Eigenvalues and Eigenvectors 



True 


Values 


Computed Values 


Eigenvalues 


Eigenvectors 


Eigenvalues 


Eigenvectors 


-11 


1 


-10.99999999 


1.0 




-29/6 




-4.833333333 




153/40 




3.82500000 


-5 


10/33 


-5.0 


! 0.3030303030 




1 




1.0 




-27/154 




-1.753246753 


-2 


0 


-2.0 


0.0 




1 




1.0 




0 




0.0 


3 


0 


3.0 


-0.0000000001 




1 




1.0 




-5/6 




-0.8333333333 


4 


1 


4.0 


1.0 




0 




0.0 




0 




0.0 


9 


0 


9.0 


0.0 




0 




0.0 




1 




1.0 



Total running time: 3.85 seconds 

Average number of iterations per eigenvalue: 4.0 
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Example 6 



This final example of the generalized problem is constructed so 
that it can also be solved as a standard problem. It is defined as follows, 
i.) A 5 (B+zC) = 0, 

where B and C are 12x12 matrices with constant elements. By multiplying 
on the left by C ^ a transformation is made to the standard problem, 



This problem has been solved in both forms by the technique discussed here. 
In addition, it was solved using the algorithm in (5), so that the results 
might be compared. 

The matrices B and C are defined on the following page. 



ii.) A' 3 (BC“VzI) = 0. 
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The matrices B and C were defined as: 



f 10 


-15 


1 


-2 


-81 


22 


-46 


-25 


-85 


-3 


24 


48 


-22 


-97 


-76 


-42 


-93 


-6 


61 


7 


37 


-39 


-81 


16 


-6 


-77 


-6 


11 


-42 


27 


99 


-72 


56 


69 


9,8 


-96 


-91 


-5 


7 


18 


-89 


14 


-63 


-10 


-17 


-85 


-36 


53 


■ 53 


53 


28 


69 


-88 


33 


-70 


^-63 


-40 


-48 


-3 


49 


1 — 1 
00 

1 


-72 


-4 


-96 


24 


29 


-20 


68 


26 


-46 


0 


57 


39 


66 


-84 


5 


-4 


-25 


26 


-44 


-91 


26 


-64 


-94 


26 


0 


-4 


-87 


-77 


42 


0 


69 


-62 


56 


86 


-69 


-65 


95 


-55 


18 


25 


57 


29 


88 


-67 


-9 


-83 


-73 


12 


-30 


-91 


-42 


27 


30 


4 


L -17 


-56 


-90 


-49 


20 



-69 


14 


62 


-36 


-20 


99 ^ 


27 


■53 


-93 


-34 


52 


19 


-15 


24 


49 


-32 


30 


-19 


3 ^ 


; -53 


-71 


57 


0 


-74 


60 ^ 


-81 


49 


60 


14 


-6 


18 


-70 


-90 


-15 


21 


-81 


71 


18 


-44 


48 


-63 


-21 


-94 


-56 


69 


-60 


-18 


-84 


57 


84 


-25 


12 


58 


-44 


-38 


62 


8 


-17 


-16 


11 


-56 


-5 


90 


-31 


-1 


85 


-18 


-72 


52 


-20 


12 


- 90 , 


66 


14 


-76 


14 


13 


00 


89 


81 


-86 


-12 


-92 


57 


32 


61 


-98 


-96 


64 


64 


-37 


63 


45 


66 


75 


66 


39 


-22 


71 


64 


91 


42 


-74 


-99 


-81 


-42 


43 


-76 


76 


36 


-84 


-93 


76 


-65 


-26 


-8 


-40 


59 


-29 


80 


18 


-82 


65 


69 


-33 


54 


28 


-35 


5 


41 


34 


37 


29 


99 


55 


-84 


29 


9 


-6 


20 


-18 


2 


73 


- 83 , 



91 

89 

-64 

16 

91 

53 

31 

20 

-18 

-59 

79 

-69 

-82 

-20 

-40 

44 

25 

35 

88 

27 

48 

-18 

-86 

-59 
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These three cases were run with a relative error bound of 10 and 
the results from all of them were identical to six or seven significant 
figures. The eigenvalues found are listed in the table below. 



Table of Eigenvalues 
.956212 
-2.67567 

- .0213507 

- .466901 

- .626829 + 2.56009 

- .274574 + .327271 
.254714 + .457145 

- .796302 + 1.04215 



Total time for present method: 104 seconds 

Average number of iterations per eigenvalue: 17.2 
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V. Areas of Further Investigation 



Some additional features are being considered for inclusion in the 
present code. However, at the time of this report all problems involving 
these revisions had not been resolved. 

It is felt that some process which could quickly determine the 
approximate location of eigenvalues may be useful, especially when large 
order problems are undertaken. It has not, as yet, been determined what 
method might be most profitable for this purpose. 

An automatic entry into double precision would be desirable for parts 
of the computations. But it will be necessary to devise some general test 
to indicate when this move could best be made without wasting computer 
time. 

A test to find when the domain of indeterminancy has been reached in 
the search for a root would be a valuable addition to those existing 
tests for acceptance of a root. 
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APPENDIX III 



A. IDENTIFICATION; 



Title: 



A Method for the Solution of the "Generalized 
Eigenvalue Problem" witti^ Polynomial Matrix Elements, 



CO-OP ID: 
Category: 
Programmer: 
Date: 



F4-NPGS-GENEIG (1604 F-63) 

Mathematical Subroutine 
R. D. Brunell 
June 10, 1966 



B. PURPOSE: 



To find those values of the complex parameter, z, which 
satisfy the system of equations: 

A^^(z)x^ + + ... + A^^(z)x^ = 0 



A2i(z)Xi + + ... + A2^(z)x^ = 0 



( 1 ) 



A T (z)x + A „(z)x + ... + A (z)x = 0, 

nl 1 nz z nn n 



where the j * ^ complex polynomials of degree 2 in z. 

This may be written more concisely as 



(2) H(z)X = 0, 

where H(z) is an n x n matrix whose elements are the polynomials 
in z, and X is a complex column vector of n elements. The elements 
of the set £z.:H(z.)X. = oj are called the eigenvalues, and the 
X^ are called the associated eigenvectors. 

The matrix H(z) , because of the form of its elements, may be 
wisitten in the following form: 

(3) H(z) = A^z^ + A^z + A^ , 

where the n x n matrices A^ have constant complex elements. 

Equation (3) will be useful later when discussing the input for the 
subroutine. 
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C. USAGE: 



1. Set up for GENEIG 

GENEIG is a subroutine and the user must supply a main 
program which will set up the input data for GENEIG* All 
communication to and from this subroutine is done through 
common storage, and the following DIMENSION, COMMON and TYPE 
statements must be declared by the user. 

DIMENSION A(20,20,3) ,EIGG(40,2) ,EVAL(40) ,EVEC(20,40) 

COMMON A, EIGG , EVAL , EVEC , NN , lORD , IE , 

NEIG , IMAX , ISAC , IPRINT , IVEC , EPI , EP2 , EP3 

TYPE COMPLEX A, EIGG, EVAL, EVEC 

2. Input for GENEIG 

The three-dimensional array, A, must contain the elements 
of the matrix H in the manner described below: 



a°[. = A(I,J,1), a^[. = A(I,J,2) and . = A(I,J,3) , 
ij 



(k) 

where a . . is the ij-th element of A, . 

ij k 



The two-dimensional array, EIGG, must contain the estimated 
valile for the eigenvalues as follows: 

EIGG(I,1) will be considered as the estimate to the value of the 
i-th eigenvalue, and EIGG(I,2) will be considered as the probable 
precentage error in that guess. The first three values necessary 
to start the process for the i-th value are computed by 

= EIGG(I,1) - EIGG(I,2)*EIGG(I,1) 



(4) 



( 2 ) 

z ^ = EIGG (1,1) 

(3) 

z ^ = EIGG(I,1) + EIGG(I,2)*EIGG(I,1) 



The array EVEC will contain the associated eigenvectors. The 
eigenvector associated with the i-th eigenvalue will be found in 
EVEC (J, I), where J = 1,2,. ..,NN. 
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NN is the order of the coefficient matrices. 



lORD is the degree of the polynomial elements. 

IE is the number of the first eigenvalue to be found. 

NEIG is the number of the last eigenvalue to be found. 

IMAX is the maximum number of iterations to be allowed. 

ISAC is a flag for initial scaling. If ISAC = 0 no scaling 
will occur. 

IPRINT is a flag for printing at every iteration. If IPRINT = 0 
only the final results will be printed. 

IVEC is a flag for the eigenvector option. If IVEC = 0 
no eigenvectors will be computed. 

EPI is the relative error bound on the eigenvalues. If 

|z^ - z^_^| / |z^| 4 EPI the problem is considered to have converged. 

-4 -8 

A value of 10 ^ EPI ^ 10 has been found to be satisfactory 

for most cases. 

EP2 is an absolute error bound on the function. If |h(z^)/ ^ EP2 

the problem is considered to have converged. A value of EP2 10 

appears to be satisfactory for those cases for which this test becomes 
effective. 

EP3 is a relative error bound on the eigenvector. If B is the maximum 

element of H(z^), then if | H(z^)xJ ^ |EP3 *b| the eigenvector 

-4 

is considered acceptable. A value of EP3 ^10 is satisfactory. 



3. Calling Sequence: 

Once the common storage has been set up in the user’s program 
it is only necessary to insert the instruction CALL GENEIG. The 
subroutine will then compute all eigenvalues and eigenvectors, that 
have been indicated by the input, before returning control to the 
calling program. 
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4. Space Required: \q excluding common storage 

5. Common Storage Required; 4251^^ 

6. Print-Outs: If IPRINT = 0 then at each iteration the number of 

the iteration, the present value of and H(z^) will be printed. 

When an eigenvalue is accepted, a message to this effect will be 

printed along with the accepted value and the function value. 

If IVEC = 0 and an eigenvector is found it will be printed. 

However, if no acceptable eigenvector is found, a statement to 

this effect will be printed along with the value of the final attempt. 

7. Error Stops: If lORD is found to be other than an integer 1 or 2 

the program will print the value of lORD and stop (Exit to Monitor) . 

If the determinant evaluator is unable to complete its operation the 

code will indicate that this has happened and stop (Exit to Monitor) . 

8. Input and Output Formats: Not applicable. 

9. Input and Output Tape Mountings: Output is to Logical unit 51 

(Standard on CO-OP Monitor). 

10. Selective Jump and Stop Settings: Not applicable. 

11. Timing Examples: 



Case 


1 


2 


3 


, 4 


5 


. 6 


Size of Matrix 


3x3 


4x4 


5x5 


8x8 


12x12 


16x16 


No. of Eigenvalues 


6 


4 


5 


8 


12 


4 


No. of Eigenvectors 


6 


4 


0 


0 


0 


0 


Average No. of Iterations 


4.0 


8.25 


6.0 


10.5 


18.7 


17.5 


Total Time (seconds) 


3.85 


3.85 


4.38 


16.0 


104. 

, 


55.28 



12. Accuracy: 

Results from this code have been compared with results taken from 
several papers on the standard eigenvalue problem and they have been 
as accurate or more accurate on all tests. Also, several problems 
were run for which exact results were known. On problems with well- 
defined roots nine or ten figure accuracy was observed. When pathological 
roots were encountered the accuracy of this code compared very favorably 
with other methods tested. 
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13. Cautions to User: The DIMENSION, COMMON and TYPE statements 

as shown in Section C must be used exactly as indicated. 

14. Equipment Configuration; Standard CO-OP Computer with FORTRAN 63 
Compiler and CO-OP Monitor System. 

15. References: 

(a) R. D. Brunell, An Iterative Solution to the Generalized 

Eigenvalue-Eigenvector Problem ^ Technical Report/Research 
Paper No. 69, U..S. Naval Postgraduate School, Monterey, 
California, 1966. 

(b) J. F. Traub, Iterative Methods for the Solution of Equations , 

Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1964. 

(c) I. Tarnove, ’^Determination of Eigenvalues of Matrices Having 

Polynomial Elements,” Journal of SIAM , June 1958. 



D. METHOD : 

Only a general discussion of the method will ‘be presented here. 
For a more detailed description see Reference (a). 

The heart of the method used in this subroutine is the iterative 
root-finding technique, which was developed by J. F. Traub (b) . In 
essence the technique is to perform a quadratic fit using the last 
three iterants and their functional values. The new estimate to the 
root is defined as that root of the fitting quadratic which is closer 
to the last iterant. The process converges with almost arbitrary 
estimates; however, through the use of good estimates a considerable 
saving in computer time can be realized. Fewer iterations will be 
necessary and only those eigenvalues of interest need be computed. 

In order to avoid convergence to a previously-calculated root 
a method of suppressing zeros, first suggested by G. E. Forsythe 
(Reference (c) page 164), was incorporated into the code. This 
feature can also be used to eliminate re-calculation of roots found 
on a previous run, especially when they may lie pathologically close 
to a desired root. This can be accomplished by loading the m known 
roots into the first m positions in the EVAL array and setting IE 
equal to m + 1. 

To alleviate the problem of floating-point exponential overflow 
a two-part, a priori and a posteriori ^ automatic scaling procedure 
was devised for this code. 
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The a. priori scaling is optional and is controlled by the 
flag ISAC. This scaling occurs only once and is done prior to 
any calculations for eigenvalues. The coefficients of each 
polynomial element are summed and a maximum ,B, is found in each 
row. The scaling is carried out by dividing each coefficient of 

every polynomial in a row by where 2 ^ 

The a. posteriori scaling is not optional and will occur only 
when an overflow is sensed. At this point the problem is trans- 
formed so that the scaled eigenvalue will be exactly the original 

k+1 k k+1 

eigenvalue divided by 2 , where 2 z. 2 and z. is the value 

of the iterant at the time of the overflow. 

The computation of the associated eigenvectors is an optional 
feature in this code. The basic assumption in this section is that 
the eigenvalue reduces the rank of the original matrix by 1. Under 
this assumption 2NN attempts will be made to find an acceptable 
eigenvector. If all attempts fail, a printout to this effect will 
be made. 
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